Living on the edge of chaos: minimally nonlinear 

O models of genetic regulatory dynamics 

o 

CN Rudolf Hanel^, Manfred Pochacker^ and Stefan Thurner^'^ 

^^ ^Section for Science of Complex Systems, Medical University of Vienna, Spitalgasse 23, 

, ^ A-1090 Vienna 



l> 



C/2 



a 



I 

o 
o 



in 

o 
o 



>$ 



^Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA 
E-mail: thurnerOunivie . ac . at 



(-H Abstract. Linearized catalytic reaction equations - modeling e.g. the dynamics of genetic 

(~| regulatory networks - under the constraint that expression levels, i.e. molecular concentrations 

of nucleic material are positive, exhibit nontrivial dynamical properties, which depend on the 
average connectivity of the reaction network. In these systems the inflation of the edge of chaos 
and multi-stability have been demonstrated to exist. The positivity constraint introduces a 
nonlinearity which makes chaotic dynamics possible. Despite the simplicity of such minimally 
nonlinear systems, their basic properties allow to understand fundamental dynamical properties 
of complex biological reaction networks. We analyze the Lyapunov spectrum, determine the 
probability to find stationary oscillating solutions, demonstrate the effect of the nonlinearity on 
the effective in- and out-degree of the active interaction network and study how the frequency 
distributions of oscillatory modes of such system depend on the average connectivity. 



^_^ 1. Introduction 

Q\ Many complex systems in general - and living systems and cells in particular - display remarkable 

^~~' stability, i.e. a capacity to sustain their spatial and temporal molecular organization. Yet, their 

stability is dynamic, i.e. these systems - to a certain degree - are capable of adapting to 
changes in their physical and chemical environment. This has led several authors [1, 2, 3, 4] 
to interpret such systems as existing at the edge of chaos. Mathematically the edge of chaos 
refers to regions in parameter space, where the system dynamics is characterized by a maximal 
Lyapunov exponent (MLE), Ai, equal to zero. In this case small changes in parameters may 
cause the dynamics to switch between regular and chaotic behavior, thereby being able to explore 
large portions of the system's phase-space. This possibility is most relevant for living systems 
C^ existing in fluctuating environments. In many dynamical systems the edge of chaos exists only 

for a tiny portion of parameter space, typically in sets of singular points, i.e. sets of measure 
zero. The dynamics of systems at the edge of chaos can become highly nontrivial, even for 
simple maps like the logistic map [5]. It has been argued that living systems have evolved 
towards the edge of chaos by natural selection [2], however it is not clear which mechanisms 
allow self-organization around these exceptional regions in parameter space. 

Living systems have exist in the state of quasi-stationary nonequilibrium and therefore can 
not be closed systems. They require a flow of substrate and energy to and from the system. 
Since long [6], rate equations for molecular dynamics have been considered. For systems to be 
self-sustaining, such rate equations need to be autocatalytic, i.e. some molecular species directly 



or indirectly catalyze their own production. For living systems, cells in particular, to be in a 
stationary state, production, decay and flow rates of intercellular components effectively have 
to balance each other, [7, 8]. Replicating, living systems therefore in general balance between 
stationary states (nonreplicating modes) and growth (replicating mode), limited by constraints 
posed by the environment. This balance provides a natural selection criterion. 

Autocatalytic systems are frequently governed by nonlinear equations for enzyme-kinetics, 
e.g. Michaelis-Menten differential equations [9], or more general replicator equations, see e.g. 
[10]. For various reasons linearized autocatalytic networks have been considered, for the case of 
abundant substrate, see e.g. [11, 12], or for reverse engineering [13, 14]. Systems with linearized 
dynamics can be easily depicted in terms of directed reaction networks, where nodes represent 
molecular species. Two nodes, where one node directly influences (production or inhibition) 
the other, are connected by a directed link. Weights of such links quantify associated reaction 
rates; negative rates indicate inhibitory links. Weights of self-loops in the reaction network, i.e. 
links of a node onto itself, quantify decay rates. Recent progress in genomic and proteomic 
technology begins to reveal facts about regulatory networks of organisms. There is some 
evidence that these directed networks show scale- free topological organization [15, 16, 17, 18]. 
More recent evidence suggests topological differences between in- and out-degree distributions 
[19, 20]. Basically two main approaches for modeling catalytic networks have been pursued: 
Discrete state approaches, e.g. Boolean networks [21], and continuous approaches, relying on 
ordinary or stochastic differential equations [22, 23, 24, 25, 26]. The relevance of noise has been 
experimentally demonstrated [27, 28, 29, 30, 31]. 

Interestingly, various models of disordered recurrent networks [21, 32] seem to share three 
distinguished modes of operation: (i) stable, (ii) critical, and (iii) chaotic super-critical. These 
properties could be generic or even universal. The importance of determining the minimum 
complexity of models exhibiting these properties has been pointed out [21] and the question 
has been raised whether these properties can already be found in linear systems. Following 
this philosophy we have recently introduced a model for genetic regulatory dynamics [33]. This 
model is governed by sets of linear equations 
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where Aij is the weighted adjacency matrix of the full autocatalytic reaction network, whose 
entries may be zero, positive and negative - indicating that i either has no influence on j or 
the production of molecular species i is stimulated or suppressed by j, respectively. This means 
that if substrate j exists, i gets produced (or reduced) at rate Aij. Xi is the concentration of the 
molecular species i (e.g. proteins or mRNA). Jj corresponds to a flow- vector. The molecular 
species i flows into the system if Jj > and out of the system if Jj < 0. Ui is a suitable noise 
term. Negative molecular concentration values Xj do not make any sense, hence we impose the 
positivity condition 

Xi>0, for aU i . (2) 

In particular if Xj = and Eq. (1) gives Xi < 0, then effectively Xi = 0. Therefore, the 
concentration Xi will remain zero until Eq. (1) gives Xi > 0. We refer to this as the minimally 
nonlinear (MNL) model. 

MNL models have nontrivial properties [33]: (i) They have a possibility for chaotic 
dynamics, (ii) MNL exhibit an inflated edge of chaos. The positivity condition causes the small 
neighborhood of a singular point in parameter space (linear system without positivity condition) , 
with MLE Ai ~ 0, to form an extended region (plateau). This effect gives random strategies of 
evolutionary phase-space sampling a finite chance of locating this particular region in parameter 
space. This may offer an explanation for why and how complex chemical reaction systems may 



have found the vicinity of the edge of chaos at ah, before evolutionary self-organization could 
take over for an eventual fine tuning, (iii) MNL models show multi- stability. Perturbations (or 
moderate noise levels) can push the system from one attractor of the dynamics to another. 

These facts raise interesting questions, (a) The existence of chaotic dynamics in MNL systems 
straight forwardly suggests to analyze the Lyapunov spectrum of the dynamics, which encodes 
information about the attractor of the dynamics. Numerical simulations indicated so far that 
MNL systems exhibit several properties, which are particularly interesting for modeling living 
systems, (b) Can topological differences between in- and out-degree distributions [19, 20], be 
explained by MNL dynamics? MNL dynamics can down-regulate the concentration Xi of a 
fraction of nodes i to zero. These nodes i then cease to play an active role in the dynamics 
of the MNL system. The remaining nodes continue to play an active role in the dynamics 
and constitute the active regulatory reaction network. The active network may have topological 
properties that differ from the full network, (c) How probable is it to find oscillating dynamics in 
MNL systems, and how are fundamental frequencies of oscillatory dynamics distributed? MNL 
dynamics frequently shows oscillatory dynamics. This is particularly interesting, since periodic 
dynamics are well known in regulatory networks in the context of the cell-cycle, e.g. [34], or 
circadien clocks [35]. Evidence has been presented that oscillating regulatory networks are also 
involved in the morphogenesis of mice [36]. Moreover, eukaryotic cells may encode information 
about extracellular environment in the frequency of stochastic intracellular events, rather than 
in the concentrations of molecular species [37] . Intracellular dynamics in terms of (stochastic) 
rhythmic burst, may be a common mechanism of intracellular information transduction . 

In Section 2 we give a summary of the model [33]. In Section 3 we report results on the 
properties of MNL systems. In Section 4 we conclude. 

2. The stochastic MNL model 

We present the MNL model as introduced in [33]. There we derived Eq. (1) by linearizing a set 
of nonlinear differential equations 

where the state vector y represents a collection of concentrations yi of molecular species i. 
These molecular species include both mRNA and proteins. The state vector y can be written as 
y = (xi, . . .Xn,Pi, ■ ■ -Pm), where the Xi, with i = 1 . . .m, are concentrations of mRNA and the 
Pr, with r = l...m, are protein concentrations. Equation (3) can be linearized around some 
fixed point y^ = (x^, . . .x^,Pi, . . .p%). The variables pr can be eliminated by the assumption 
that, around a fixed point, changes of mRNA (or protein) concentrations, 5x = x — x", 
translate linearly into variations of the protein (or mRNA) concentrations, Sp = p — p^, i.e. 
^Pr = "^iCriSxi, for some fixed matrix C. Assuming (thermal) fluctuations of the production 
and degradation rates around average values Aij, using the law of large numbers, finally leads 
to the equation 

3 

where ^j G A^(0, &) and r/j G -^(0, cr) are independent normally distributed zero-mean random 
variables with standard deviations a for the multiplicative noise and a for the additive noise. 
Comparing with Eq. (1) the noise-term vi can be identified, Vi = E,i{t){xi — x?) -|- ryj, and the 
flow-vector J is 

Ji ^ — y ^ AijXj . (5) 

j 

Time series of mRNA expression levels Xj(t) typically oscillate around fixed points (average 
values) x" = {xi{t))t. This can be used to directly feed characteristic mRNA expression profiles 



:Xi = Y, Aij {xj - x°) + ^i{t) (xi - x°) + r,i , (4) 



into the MNL model. Although from a purely mathematical point of view, the fixed point x*^ = 
is a perfectly legitimate choice, it contradicts the fact that living systems are open systems and 
require non-vanishing effective flow-vectors, J 7^ 0. Equation (5) immediately implies that the 
choice x'^ = is incompatible with this requirement. Here we use x^ = 1000, for all i. 

The weighted adjacency matrix Aij can be used to incorporate topological information on 
biological networks. In our model the decay rates. An < 0, have identical value, for all i. 
This assumption is wrong in general but reasonable for niRNA encoding groups of proteins 
that act together in stoichiometric complexes [38]. Since A is largely unknown experimentally 
we are interested in random ensembles of matrices A, which can be parametrized with only a 
few parameters. We model ^ as a random matrix in the following way. Using terminology 
from network theory, the out-degree ki of a node (= molecular species) i is defined as the 
number of products j that can be regulated by the molecular species i. The in- degree is the 
number of molecular species that regulate i. The ensemble of interaction networks can now 
be specified by the in- and out- degree distribution Pm/outi^)- Although in principle in- and 
out-degree distribution can be chosen independently, we consider identical in- and out-degree 
distributions. In [33] we have compared scale- free networks with Erdos-Renyi networks [39] 
and have noted only minor effects on the stability, i.e. the formation, of the Ai ~ plateau. 
Here we will only consider Erdos-Renyi networks. The associated topological ensembles for the 
full reaction networks are completely specified by the number A^ of molecular species and the 
number L of links between them, i.e. the average degree {k) = L/N of the network. 

Once the topology of a network is fixed, the actual weights Aij are sampled from a 
normal distribution N{0,aA) with zero mean and standard deviation aA- This assumption 
is experimentally supported by e.g. [40]. We define the constant D = —An/ a a and set cja = 1 
in all numerical simulations. The time-increment used for all numerical simulations is dt = 0.1. 

The maximal Lyapunov exponent, Ai, measures the exponential rate with which a 
perturbation of a trajectory propagates over time. If Ai < 0, the perturbation vanishes. If 
Ai > 0, the perturbation grows exponentially. In [33] it was shown that for MNL systems an 
interval of average network connectivities / = [/c~. A;"*"] exists so that Ai < for {k) < k~ (Area 
A), Ai ~ for (k) £ I (Area B), and Ai > for (k) > /c"*" (Area C). The two values of (k), which 
estimate the beginning and the end of the Xi{{k)) ~ plateau, i.e. the interval /, are given by 
k~ = D^ and k^ = 2D^. It also was shown, that as (A;) gets larger than D^, the number A^zcro 
of nodes i, whose concentration Xj = increases monotonically, due to the positivity condition, 
until A^zero ~ N/2. The sub-network of size A^on = N — N^ero, consisting of those nodes j of the 
full network, which have nonzero concentration Xj, we call the active network. If two nodes, i 
and J, of the active network have a link Aij in the full network, then the active network inherits 
this active link. We denote the adjacency matrix of the active network with A°'/. 

Technically, the positivity condition is implemented so that Xi{t+dt) = 0, \ixi{t)+Xi{t)dt < 0. 
In [33] the positivity condition, Eq. (2), was implemented in a slightly different way, i.e. 
Xi{t + dt) = Xi{t), if Xi{t) + Xi{t)dt < 0. The different implementation has a small effect 
on the plateau formation. 

3. Results 

We now present (i) the Lyapunov spectrum of the MNL model, (ii) the probability to find 
growing, decaying, or stable dynamics and the size and topological properties of the active 
catalytic network. Further we present (iv) the probability of finding oscillating time series and 
their characteristic frequencies. In all following figures the theoretical plateau interval [k~ , k'^] 
is marked (gray shading). 
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Figure 1. (a) Largest ten Lyapunov exponents (Ap, p = 1, . . . , 10) of the Lyapunov spectrum 
(A^ = 500). The inset magnifies the plateau region. The two black dashed lines are theoretical 
curves, derived in [33], approximating Ai((/c)) in the areas A and C. The intersection of these 
curves with the x-axis, Xi{{k)) = 0, estimate the beginning and end of the Ai((fc)) ~ plateau 
(area B). (b) The Lyapunov spectrum as a function of (k) is shown (A^ = 200). (c,d,e) Lyapunov 
spectra, \p{{k)) as functions of p, for areas A, B, and C. Li area A spectra are all similar, in 
area B the slope of Xp gets steeper with growing (k), while Ai((/c)) ~ 0. Comparing (d) and 
(e) shows a clear qualitative difference between area B and C. The Kolmogorov Sinai Entropy 
EkSj diagram (f), and the Kaplan York Dimension Dky, diagram (g), both divided by system 
size A^, in logarithmic scale, for N = 200 and A^ = 500. Dky represents an upper bound for the 
information dimension of the system. Eks is calculated from Pesin's theorem, as the sum over 
all positive Ap > in the Lyapunov spectrum, for fixed (k). 

In (a-g) averages are taken over 50 random realizations, time interval [200, 1000], noise 
a = a = 0.1, CJA = 1, £> = 4, and x° = 1000. 



3.1. The Lyapunov spectrum 

It is straight forward to compute the fuh Lyapunov spectrum Ap which ahows to determine 
properties of attractors in greater detail. In particular we computed the Kaplan- York Dimension 
Dry which gives an upper bound for the information dimension of the system and the 
Kolmogorov- Sinai Entropy Ers^ see e.g. [41]. While Dry gives an estimate of the dimension 
of the attractor, i.e. the phase-space volume-preserving subspace of the dynamics, Ers can be 
interpreted as a measure of the number of excited states in the system. 

In Fig. (1) we summarize numerical results, (a) shows the first ten Lyapunov exponents 
Xp, p = 1 ... 10, as functions of (k). Clearly, Ai ~ matches the theoretical plateau region, 
[k~ , k~^]. (b) shows how the full Lyapunov spectrum depends on (k). In area A {{k) < k~) all 
Ap are densely arranged. In area B (the plateau) Ai ~ 0, while all Xp decrease monotonically for 
p > 1. As a consequence, the Lyapunov spectrum gets less dense with growing (k) and covers 
an increasing range of negative values. In area C {{k) > fc"*") the Lyapunov spectrum gets still 
less dense. Yet, this decrease in density is qualitatively different than in area C. Ai is increasing 
in area C, while Ap, for large p, still decreases monotonically, but less pronounced than in area 
B. This can be seen clearly in Fig. (1) (c-e), where Lyapunov spectra for various values of (k) 
are shown as functions oi p. In (c) the values of (k) are chosen from area A, in (d) from area B, 
and in (e) from area C. 

Further, Fig. (1) shows the average Kolmogorov-Sinai entropy (f) and the average Kaplan- 
York dimension (g). Both quantities require the existence of some Ap > 0, i.e. Ers and 
Dry cannot be computed for area A. In area B and C averages of Ers and Dry can only 
be taken over realizations, which have at least some Ap > 0. Clearly, in area C the entropy 
per node Ers/^ and the fraction Dry /N of the volume-preserving subspace seem to become 
independent of system-size N , for sufficiently large (A;) > k'^ . However, in the plateau region, 
area B, both quantities do not seem to scale with system size and finite size effects may become 
relevant. It is an interesting open question towards which limit-function Eks/^ and Dry /N 
converge as A^ — )• oo. 



3.2. Stability of MNL systems 



Stability 




Figure 2. Fraction of realizations which lead to exponentially growing (Ai > 0.1), decaying 
(Ai < —0.1) and stable time series (|Ai| < 0.1) computed from 100 realizations, N = 500, D = 4, 
time interval, [200, 1000], a = a = 0.1, aA = 1, and x°=1000. 

What is the probability of finding the dynamics of MNL systems to be characterized (i) 



by exponential growth, (ii) exponential decay, or (iii) non-exponentially growing stationary or 
oscillatory dynamics? Living systems can be expected to exist close to stationary or oscillatory 
states [7, 8]. Sufficiently positive and sufficiently negative MLEs, Ai, in complete analogy to 
linear systems, indicate exponential growth or decay. Therefore the probabilities of finding MNL 
systems in one of the growth modes (i-iii) can simply be estimated by thresholding Ai, for a 
sufficiently small threshold e > 0. Counting realizations in the MNL ensemble with (i) Ai > e, 
(ii) —E > Ai, and (iii) e > |Ai| estimates the ratios of growth mode fractions (i-iii)- Numerical 
results for e = 0.1 are given in Fig. (2). Clearly, dynamics of type (iii) is favored in the plateau 
region. Yet, to a much lesser extent, stationary and oscillatory dynamics also can be found in 
areas A and C. 

3.3. The active genetic regulatory sub-network 
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Figure 3. Average fractions of Xj: positive 

are taken over 1000 realizations, time interval, [500, 1000], A^ 

xO = 1000. 
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Recent evidence from the analysis of genetic regulatory networks [19, 20], suggests topological 
differences between in- and out-degree distributions. Can differences between in- and out-degree 
distributions appear merely by the fact that the full interaction network Aij is different from 
the active sub-network ^°?? Is it possible that through a symmetry-breaking mechanism the 
active in- and out-degree distribution pf^{k) and Pout(^) become different from p^'^^^(/c)? 

We have analyzed the average properties of active sub-networks of the MNL model and 
distinguish three types of nodes in MNL systems: (i) nodes i with concentrations Xj > for all 
times, (ii) nodes j with Xj = for all times, and (iii) nodes that alternate between on and off. 

zero + A^ait and 



The associated numbers of nodes are Aposi Azero, and A'ait, where A^ 



the fraction of nodes are denoted n 
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Npos/N, 



iVpos + N^, 



riy. 



N^cro/N, and rzait = N^it/N. Note 



A^pos + A^ait- In Fig- (3) we show these fractions 



that the size of the active sub- network is Aon 
for a network with A^ = 500. Alternating nodes are most important in the plateau region, where 
A^zero starts to grow, i.e. the active networks shrinks with growing (k). For (k) > k^ the fraction 
Ti-ait decreases and reaches a constant value riait ~ 0.05; ripos and rizero become equally large. 

In Fig. (4 (a) unweighted in- & out-degree shows the in- and out-degree distributions of the 
active network for various {k) > k~ . The degree distribution of the full network is shown (red) 
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Figure 4. (a) Unweighted in- and out-degree distributions of the active regulatory sub-network 
for various (A;). Active in- and out-degree, P°n/out(^)) are practically indistinguishable, (b) 
Weighted in- and out-degree distributions. In- and out-weight distributions, P°°/out(0)> of active 
weights are clearly distinguishable. 4> = Yl ^ij ^^^ ^^^ ^^^^ runs over i or j for in- and out-weight 
distribution, respectively, (c) Mean, (d) standard deviation, (e) skewness and (f ) kurtosis of the 
in/out-weight distributions. Differences between in- and out-weight distributions are found in 
the standard deviation and the skewness. Averages are taken over 50 realizations, N = 500, 
time interval, [500, 1300], D = 4, a = a = 0.1, a^ = 1, x° = 1000. 



for reference. Although in- and out-degree distribution of the active network differ substantially 
from the degree distribution of the full network, in- and out-degree distributions essentially 
remain identical. If we look at the weight distributions, p°^ and PoUt' associated with active in- 
and out-links in (b) weighted in- & out-degree the situation changes: differences in the in- and 
out-weight distributions begin to show. These differences are recognizable in (d) the standard 
deviation and (e) the skewness of the weight distributions, but not in (c) the mean and (f ) the 
kurtosis of the active weight distributions. This establishes evidence that a possible symmetry 
breaking of in- and out-degree distributions of complete chemical reaction networks can arise 
due to the natural nonlinearity in the dynamics of the chemical reactive systems, i.e. Xj > for 
all i, at all times. However, the size of this effect seems to be insufficient to explain the size of 
topological differences [19, 20]. 

3.4-- Oscillating modes in MNL systems 
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Figure 5. (a) Probabil- 
ity of finding oscillating 
realizations: with ex- 
isting fundamental fre- 
quency cj* (blue) and 
both, existing ojI and 
^2 (green). (b) Av- 
erage ujI as a func- 
tion of {k). (c) Stan- 
dard deviation of wjj". 
A^ = 500, time inter- 
val, [1000, 3000], D = 
4:, a = a = 0, aA = 
1, x° = 1000. In (b) 
and (c) N = 500, D = 
6, (green circles) and 
N = 200, D = i (red 
squares) are shown for 
comparison. 
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What is the fraction of MNL systems, which display oscillating dynamics and what are their 
typical frequency distributions? We find that if a particular realization of an MNL system shows 
oscillating dynamics, then all Xj in the active network of the particular realization follow the same 



fundamental oscillation-pattern. The dominant frequencies w*, s = 1,2, 



correspond to 



local maxima in the power-spectra of the active Xj. Smax < -/V is the maximal number of 
detectable local maxima in the power-spectrum of the MNL system dynamics. We looked for 
fundamental frequencies w J (if existing) . 

Technically we identified oj^ and 0J2 in the following way. We computed the power-spectrum 
Pi{uj) for each time series Xi{t) in a particular realization. We took the weighted average 
g(u}) = {wiPi{uj))i over all nodes of the realization. The weights, Wi = Pi{uJi)~^, have been 
chosen inverse-proportional to the power Pi{u}i) of the lowest frequency oji in the power- 
spectrum. This choice turned out to be optimal for correctly detecting the dominant frequencies 
oj* of MNL dynamics. The first frequency ojI can be found in the following way. We have searched 
for the local minimum of g{uj) with the smallest frequency tDi. If no such local minimum exists 
the time series was classified as non-oscillating. If wi exists, the fundamental frequency, wf > cDi, 
is determined such that g{ujl) is the maximum of all g{u}), with uj > oJi. Similarly, we computed 
a second dominant frequency by searching for the next local minimum uj2 > 1^1, and take (7(^2) 
to be the maximum of all g{oj), with uj > u}2- If 1^2 exists, 0^2 is the second characteristic 
frequency of the system. 

In Fig. (5) (a) shows the probability of finding a realization of the MNL model possessing a 
fundamental frequency and a second dominant frequency. Oscillating realizations are dominant 
in the plateau region and the probability of finding oscillating realizations is close to certainty 
for k^ < {k) < {k^ + k'^)/2. For {k) > {k~ + k~^)/2 this high probability decreases, but still 
has a value of about 0.6 for {k) = k^ . Furthermore, in Fig. (5) (b) the average and (c) the 
standard-deviation of the fundamental frequencies ujI are shown. 

4. Conclusions 

We presented results on properties of the MNL model. We analyzed the Lyapunov spectrum 
of the model and computed the Kolmogorov- Sinai Entropy and the Kaplan- York Dimension, 
characterizing the attractors of MNL dynamics. We analyzed stability properties by computing 
the probabilities for finding exponentially growing, decaying and non-exponentially growing 
(stable) dynamics and found that stable dynamics plays a dominant role in the plateau interval, 
[k~ , k~^]. We determined characteristic fractions of concentration levels, which are always down- 
regulated to zero, are always positive, or are alternating, i.e. oscillating between zero and non- 
zero concentration levels. Nodes with alternating concentration levels are dominating in the 
plateau interval. We analyzed topological properties of the active regulatory network, consisting 
only of molecular species (nodes) with nonzero concentration levels in a given time period. We 
found no symmetry-breaking in the in- and out-degree distributions of the active regulatory 
network with respect to the full network Aij. However, we found symmetry- breaking in the 
in- and out- weight distributions of active networks. This indicates that in chemically reactive 
systems the natural nonlinearity introduced by the positivity condition, i.e. concentrations 
of molecular species can never be negative, suffices to implement a symmetry- breaking in the 
topology of the system, which can actually be measured. One may speculate if the pronounced 
differences of in- and out-degree distributions as found in living organisms, have their origin in 
symmetry-breaking mechanisms, which later could become amplified by selective evolutionary 
processes. Finally, we determined probabilities of finding oscillating dynamics in MNL systems 
and analyzed fundamental properties of their dominant frequencies. We found that oscillatory 
dynamics is most likely, in fact almost certain, for average connectivities of networks chosen from 
the plateau interval. This corresponds well to the observation that regulatory networks of living 
organisms, cells in particular, frequently show sub-networks with oscillatory dynamics. The 
properties analyzed indicate that near the edge of chaos MNL system - despite the simplicity 
of the MNL model - display many important characteristic properties, which are expected from 
living matter. 
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